
	Program Shoreline
	Dimension H(100,100),Hupld(100),JSHORE(100),Yoff(100),
	. SXser(1,100),SYser(1,100),Hseries(1,100,100),YSHORE(100),
	. Hold(100)
	open(unit=86,file='log.txt')
	ITLIMIT=150
	ITIME = 0
	I=1
	NN=0
	NNSAVE=5
	Yoff(I)=50.
	DeltaX = 500.
	DeltaY = 500.
	DT = 1.
	HHinc = 0.25
	Coef=0.3
	IMAX=1
	JMAX=12
c	H(1,1) = -8.5
c	H(1,2) = -8.1 !-7.5
c	H(1,3) = -7.6 !-7.1
c	H(1,4) = -7.2 !-6.8
c	H(1,5) = -6.5 !-6.3
c	H(1,6) = -4.0 !-6.
c	H(1,7) = -0.5 !-5.5
c	H(1,8) =  2. !-5.
c	H(1,9) =  2. !-4.2
c	H(1,10) = 2.
c	H(1,11) = 2.
c	H(1,12) = 2
	H(1,1) = -7.8
	H(1,2) = -7.5
	H(1,3) = -7.1
	H(1,4) = -6.8
	H(1,5) = -6.3
	H(1,6) = -6.
	H(1,7) = -5.5
	H(1,8) = -5.
	H(1,9) = - 4.2
	H(1,10) = -.5
	H(1,11) = 2.
	H(1,12) = 2
c
	JSHORE(I) = 10
 	Hupld(I) = 2.   !Remember this will need to be an input array
c
	Htop = 2.
	Hbot = -3.
c
c    End of the inputs
c
	DeltaH = Htop-Hbot
	tanALPHA = DeltaH/DeltaY
	TA = tanALPHA/2.
	Hav = Htop -  (DeltaH/2.)
	TA = tanALPHA/2.
	HAA = ((-Hbot)**2)/(2.*DeltaH)
	HBB =  (Htop**2)/(2.*DeltaH)
	HA = Hbot + HAA
	HB = Htop - HBB
c
	Ylamax = ABS(Hbot/tanALPHA)
	Ylbmax = ABS(Htop/tanALPHA)
c
	FMa = ABS((Ylamax)/DeltaY)
	FMb = ABS((Ylbmax)/DeltaY)
c
	Hold(I) = H(I,JSHORE(I)) ! Remember this need work to make subroutine
c
	IISAVE=0
c
c	Output the initial profle
c
	NN=1
	YSHORE(1)=((JSHORE(1)*DeltaY)+(DeltaY/2)) + Yoff(I)
	DO 3 I=1,IMAX
	SYser(I,NN)=YSHORE(1)
	SXser(I,NN)=((float(I))*DELTAX)-(DELTAX/2.)
     	DO 3 J=1,JMAX
	Hseries(I,J,NN)=H(I,J)
    3 CONTINUE
c
c     Begin the main time loop
c
	DO 800 IT=1,ITLIMIT
	IISAVE = IISAVE + 1
	I=1
	Hinc = HHinc
c
c	Adjust the upper shoreface cell (i.e. JSHORE(I)-1) so that it comes up
c	 Hbot as the shoreline advances or drops to the adjacent cell heights as
c      the shoreline retreats.  The concept here is that the response time of 
c      surf zone is so fast that it can be thought of as a fixed profile 
c      surrounded by an envelope of maximum and minimum shapes. Next cell
c	 seaward is the upper shoreface which tracks the surf zone in a damped
c	 manner.
c
c	 The movement of this upper shoreface cell depends first on whether the
c      surf zone is advancing or retreating landward.
c
	If(Hinc.GE.0.)Goto 14 ! advancing case transfers
c	  eroding case falls through

	Htarget = (Hbot + (H(I,(JSHORE(I)-2))))/2.
	Halter = Coef*(Htarget-(H(I,(JSHORE(I)-1))))
c	write(*,*)'H(I,(JSHORE(I)-1))=',H(I,(JSHORE(I)-1))
c	write(*,*)'Htarget=',Htarget,'  Halter =',Halter
	IF(Halter.GT.Hinc)GOTO12 !If Halter is more (neg) than Hinc,limit it
	Halter = Hinc
c
   12	Hinc = Hinc - Halter
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1)) + Halter
	GOTO18
c     The following is the advancing case
   14	Halter = Coef*(Hbot - H(I,(JSHORE(I)-1)))
	IF(Halter.LT.Hinc)GOTO16 !If Halter is more than Hinc,limit it
	Hchange = Hinc
  16	Hinc = Hinc - Halter
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1)) + Halter
c
c  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c	write(*,*)' H(I,(JSHORE(I)-1))=',H(I,(JSHORE(I)-1))
c  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c
   18 Continue !If(Hinc.EQ.0.)GOTO800
c
c     Now that the upper shoreface cell is adjusted the surf zone changes
c	are calculated. 
c
c     When converting to the subroutine we will need to compute a Hinc from
c     the difference of Hold(I) and H(I,JSHORE(I)). At this time we set this 
c     difference as an input and make the following statement
c
	 
	 H(I,JSHORE(I)) = Hold(I) + Hinc
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	Write(*,*)'Hi,j=',H(I,JSHORE(I)),'  HA=',HA,'  HB=',HB
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c
c
c  
c!!!!!!!!!!!!!!!!!!!	 Account for possible bluff erosion or overwashing
c!!!!!!!!!!!!!!!!!!!	  BSVol is the volume from bluff erosion; VVol is overwash
c
c!!!!!!!!!!!!!!!!!!!	BSvol = 0.
c!!!!!!!!!!!!!!!!!!!	IF(Hupld(I).EQ.Htop)GOTO 8
c!!!!!!!!!!!!!!!!!!!	IF(Hupld(I).GT.Htop)GOTO 7
c!!!!!!!!!!!!!!!!!!!	VVol = Vol2
c!!!!!!!!!!!!!!!!!!!	GOTO 8
c!!!!!!!!!!!!!!!!!	IF(YPoff1.LE.((DELTAY/2.)-Yref))GOTO8
c!!!!!!!!!!!!    7	BSVol = (Hupld(I)-Htop)*(YPoff2-YPoff1)
c!!!!!!!!!!!!!!!    8 Continue
c
c     H(I,(JSHORE(I)))
c	
c    Test whether the shoreline is advancing seaward or retreating landward
	IF(H(I,(JSHORE(I))).LT.Hold(I))GOTO24
c
c        The advancing case falls through
c
	IF((H(I,(JSHORE(I))).GT.Hav).AND.((Hold(I)).GE.Hav))GOTO23
	IF((H(I,(JSHORE(I))).GE.Hav).AND.((Hold(I)).LE.Hav))GOTO22
c
c         The advancing <Hav,<Hav case falls through
c
	Ylbold = Ylbmax*((Hold(I)-Hav)/(HA-Hav))
	Ylbnew = Ylbmax*((H(I,(JSHORE(I)))-Hav)/(HA-Hav))
c
	VolBold = (0.5)*((Ylbold**2)*tanAlPHA)	   
	VolBnew = (0.5)*((Ylbnew**2)*tanALPHA)	   
c
	VolDiff = VolBold - VolBnew
	hDiff1 = FMb*((Ylbold/Ylbmax))*(VolDiff/DeltaY)	  !!!old for new & max)**2
	hDiff2 = FMb*((VolDiff/DeltaY)-hDiff1)
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
	Write(*,*) 'hDiff1 =',hDiff1
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - hDiff2
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1)) + hDiff1
	GOTO40
c
   22	Continue ! This is the advancing Straddle Hav case
c
	Ylbold = -(Ylbmax*((Hold(I)-Hav)/(Hav-HA)))
	Ylanew = - (Ylamax*((H(I,(JSHORE(I)))-Hav)/(HB-Hav)))
c
	VolBold = (0.5)*((Ylbold**2)*tanAlPHA)
	VolAnew = (0.5)*((Ylanew**2)*tanALPHA)
c
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1))+(VolBold/DeltaY)
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - ((VolBold+VolAnew)/DeltaY)
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	+ (VolAnew/DeltaY)
	GOTO40
c
   23	 Continue ! This is the advancing >Hav,>Hav case
c
	Ylaold = - (Ylamax*((Hold(I)-Hav)/(HB-Hav)))
 	Ylanew = - (Ylamax*((H(I,(JSHORE(I))))-Hav)/(HB-Hav))
c
	VolAold = (0.5)*((Ylaold**2)*tanALPHA)
 	VolAnew = (0.5)*((Ylanew**2)*tanALPHA)
c
	diffVol = VolAnew-VolAold
	hDiff1 = ABS(FMa*((Ylanew/Ylamax))*(diffVol/DeltaY))		!!!here
	hDiff2 = ABS(FMa*((VolDiff/DeltaY)-hDiff1))
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - hDiff2
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	+ hDiff1
	GOTO40
c
c
c      There follows the retreating cases
c
   24	IF((H(I,(JSHORE(I))).GT.Hav).AND.((Hold(I)).GE.Hav))GOTO27
 	IF((H(I,(JSHORE(I))).LE.Hav).AND.((Hold(I)).GE.Hav))GOTO26
c
c	 The retreating <Hav,<Hav case falls through
c
	Ylbold = Ylbmax*((Hold(I)-Hav)/(HA-Hav))
 	Ylbnew = Ylbmax*((H(I,(JSHORE(I)))-Hav)/(HA-Hav))
c
	VolBold = (0.5)*((Ylbold**2)*tanAlPHA)
	VolBnew = (0.5)*((Ylbnew**2)*tanALPHA)
c
	VolDiff =  VolBnew - VolBold 
	hDiff1 = FMb*((Ylbold/Ylbmax))*(VolDiff/DeltaY)	!!!old for new	!!!here
	hDiff2 = FMb*((VolDiff/DeltaY) - hDiff1)
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) + hDiff2
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1)) - hDiff1
	GOTO40
c
c
   26	Continue ! This is the retreating Straddle Hav case
c
	Ylaold =(Ylbmax*((Hold(I)-Hav)/(HB-Hav)))
	Ylbnew =(Ylamax*(H(I,(JSHORE(I)))-Hav)/(HA-Hav))
c
	VolAold = (0.5)*((Ylaold**2)*tanAlPHA)
	VolBnew = (0.5)*((Ylbnew**2)*tanALPHA)
c
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1))- (VolBnew/DeltaY)
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) + ((VolAold+VolBnew)/DeltaY)
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	- (VolAold/DeltaY)
	GOTO40
c
   27	 Continue ! This is the retreating >Hav,>Hav case
c
	Ylaold = - (Ylamax*((Hold(I)-Hav)/(HB-Hav)))
 	Ylanew = - (Ylamax*((H(I,(JSHORE(I))))-Hav)/(HB-Hav))
c
	VolAold = (0.5)*((Ylaold**2)*tanALPHA)
 	VolAnew = (0.5)*((Ylanew**2)*tanALPHA)
c
	diffVol = VolAold - VolAnew
	hDiff1 = ABS(FMa*((Ylaold/Ylamax))*(diffVol/DeltaY))	!!!here 
	hDiff2 = ABS(FMa*((VolDiff/DeltaY)-hDiff1))
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) + hDiff2
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	- hDiff1
	GOTO40
c
   40	Continue
c
c	Trap the overfilling condition
c
	IF(H(I,(JSHORE(I))).LE.HB)GOTO42
	Delta = H(I,(JSHORE(I)))-HB
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1)) + Delta
	H(I,(JSHORE(I))) = HB
c	
c     The shoreline is tracked with the cell-referenced value YPshore which is
c	measured in the plus-Y direction from the right (i.e.offshore) cell face  
c
   42	YPShore = DeltaY*(1.-((H(I,(JSHORE(I)))-HA)/(HB-HA)))
c
c	 The shoreline is tracked in grid space by YOFF(I)
c
	YOFF(I) = ((JSHORE(I)*DeltaY) - (DeltaY/2.)) + YPShore
c
c      Determine if a cell shift is needed 
c
c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	write(*,*)'Hi,j=',H(I,(JSHORE(I))),'  YPShore =',YPShore
c	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	IF(YPShore.GT.0.)GOTO45
	IF(H(I,(JSHORE(I)-1)).GE.HA)GOTO44
	Delta = HA -  H(I,(JSHORE(I)-1))
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - Delta
	H(I,(JSHORE(I)-1)) = HA
	GOTO50
   44	JSHORE(I) = JSHORE(I) - 1
	GOTO50
   45	IF(YPShore.LE.DeltaY)GOTO50
	JSHORE(I) = JSHORE(I) + 1
c
c    Store the JSHORE(I) depth for the next time step
c
   50	Hold(I) = H(I,JSHORE(I))
c
c Store data from designated time steps
c
	IF(N.GT.1)GOTO 410					  !CHANGED 3/11
	GOTO 411							  !CHANGED 3/11
  410	IF(IISAVE.NE.NNSAVE)GO TO 800
  411	IISAVE=0
	NN=NN+1
	NDATA=NN
	DO 430 I=1,IMAX
	SYser(I,NN)=YSHORE(I)						  
	SXser(I,NN)=((float(I))*DELTAX)-(DELTAX/2.)
     	DO 430 J=1,JMAX
	Hseries(I,J,NN)=H(I,J)
  430 CONTINUE
c
  800 CONTINUE
c
c
	OPEN(UNIT=2,file='CONTRL.DAT')				 
	  WRITE(2,*) IMAX,JMAX,NDATA,DELTAX,DELTAY,DT
	CLOSE (UNIT=2)
c
	OPEN(UNIT=2,file='CSTMout.DAT')
	 DO 960 N=1,NDATA
	 DO 960 J=1,JMAX
  960	  WRITE(2,*)(Hseries(I,J,N),I=1,IMAX)
	 DO 965 N=1,NDATA
  965	  WRITE(2,*)(SXser(I,N),SYser(I,N),I=1,IMAX)
	CLOSE(UNIT=2)
c
c	
	end